Application of the confluent Heun functions for finding the quasinormal modes of 

nonrotating black holes 
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Although finding numerically the quasinormal modes of a nonrotating black hole is a well-studied 
question, the physics of the problem is often hidden behind complicated numerical procedures aimed 
at avoiding the direct solution of the spectral system in this case. In this article, we use the exact 
analytical solutions of the Regge- Wheeler equation and the Teukolsky radial equation, written in 
terms of confluent Heun functions. In both cases, we obtain the quasinormal modes numerically 
from spectral condition written in terms of the Heun functions. The frequencies are compared with 
ones already published by Andersson and other authors. A new method of studying the branch cuts 
in the solutions is presented - the epsilon-method. In particular, we prove that the mode n = 8 is 
not algebraically special and find its value with more than 6 firm figures of precision for the first 
time. The stability of that mode is explored using the e method, and the results show that this new 
method provides a natural way of studying the behavior of the modes around the branch cut points. 



QUASI-NORMAL MODES OF BLACK HOLES 

The study of quasinormal modes (QNMs) of a black 
hole (BH) has long history [1-7]. The reason behind this 
interest is that the QNMs offer a direct way of study- 
ing the key features of the physics of compact massive 
objects, without the complications of the full 3D gen- 
eral relativistic simulations. For example, by comparing 
the theoretically obtained gravitational QNMs with the 
frequencies of the gravitational waves, one can confirm 
or refute the nature of the central engines of many astro- 
physical objects, since those modes differ for the different 
types of objects - black holes, superspinars (naked singu- 
larities), neutron stars, black hole mimickers etc. [8-13]. 

To find the QNMs, one needs to solve the second-order 
linear differential equations describing the linearized per- 
turbations of the metric: the Rcgge- Wheeler equation 
(RWE) and the Zerilli equation for the Schwarzschild 
metric or the Teukolsky radial equation (TRE) for the 
Kerr metric and to impose the appropriate boundary con- 
ditions - the so-called black hole boundary conditions 
(waves going simultaneously into the horizon and into 
infinity) [1, 3]. Additionally, one requires a regularity con- 
dition for the angular part of the solutions. And then, 
one needs to solve a connected problem with two complex 
spectral parameters - the frequency ui and the separation 
constant E (E = 1(1 + 1) - real for a nonrotating BH, 
with I the angular momentum of the perturbation). This 
system was first solved by Chandrasekhar & Dctweiler[l] 
and Teukolsky & Press [14] and later developed through 
the method of continued fractions by Leaver [15]. For 
more recent results, see also [4-7]. 

Because of the complexity of the differential equa- 
tions, until now, those equations were solved either ap- 
proximately or numerically meeting an essential difficulty 
[1]. The indirect approaches like the continued fractions 
method have some limitations and are not directly re- 
lated with the physics of the problem. The RWE, the 



Zerilli equation and TRE, however, can be solved analyt- 
ically in terms of confluent Heun functions, as done for 
the first time in [16-19]. Imposing the boundary condi- 
tions on those solutions directly (see [13, 17]) one obtains 
a system of spectral equations (1) and (2) featuring the 
confluent Heun functions which can be solved numeri- 
cally. 

In this article, for the first time we present finding 
I and u> directly in the case for gravitational perturba- 
tion s = —2 in a Schwarzschild metric, i.e. we solve the 
RWE and TRE analytically in terms of confluent Heun 
functions and we use a newly developed method (the 
two-dimensional generalization of the Miiller method de- 
scribed in the internal technical report [20]) to solve the 
system of two transcendental equations with two complex 
variables. Then we use the epsilon method to study the 
stability of the solutions with respect to small variations 
in the phase condition. 

The results are compared with already-published ones 
and are found to coincide with at least 8 digits for the 
RWE and 6 digits for the TRE. For the first time, the so- 
called algebraically special mode n = 8 is evaluated with 
precision of more than 6 digits, and it is shown to have 
a nonzero real part. This firmly refutes the hypothetical 
relation of this mode with the algebraically special once. 
Also demonstrated is the nontrivial dependence on e of 
the first 11 modes in both cases. 



GENERAL FORM OF THE EQUATIONS 

The angular equation for both cases is the solution of 
the Teukolsky angular equation when there is no rotation 
(a = 0): 

S(6) = (cos(0) - l)(cos(0) + l)LegcndreP(Z, 2, cos(0)) = O 

(1) 

where 9 £ [0,7r] is the angle. The results for the QNMs 
should be independent of the choice of 9 in the spectral 
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conditions. In our numerical experiments, we use 8 = 
7T - 1CT 7 . 

The general form of the radial equations is obtained 
from the solutions of the RWE and TRE written in terms 
of the confluent Heun functions according to [16], on 
which the black hole boundary conditions have been im- 
posed. The choice of the local solution in terms of the 
Heun function takes into account the boundary condition 
on the horizon. Then, it remains to impose the following 
boundary condition on the space infinity (for details see 
[13, 16]): 

R = rP HeunC(a, p, 7, 6, T}, 1 - »•«,) = 0, (2) 

where HcunC is the confluent Heun function as defined in 
maple and the parameters a, (3, 7, S, rj and p differ for the 
two equations. The values of the parameters when the 
BH mass is M = 1/2 and, if we choose \roo \ = 20 which 
turns out to be large enough to simulate numerically the 
actual infinity, are ([16, 19]): 

1. for the solutions of the Regge- Wheeler Equation: 

a = —2 iu, (3 = 2 iu, 7 = 4, 

5 = -2u 2 ,r] = 4:-l-l 2 + 2u 2 , 

roo = 20 c- l{1 ' 2 (i+«)t+«bM) , p = 3, 

2. for the solutions of the Teukolsky Radial Equation: 

a = -2 iu, (3 = 2 + 2 iu, 7 = 2, 

S = -Aiu - 2u 2 , rj = 2u 2 + Aiu - A, 

r oo = 20 C -^ 2 ( l+ ^ +ar ^\p = b, 

where A = 1(1 + 1) — s(s + 1) is the separation 
constant. The parameters were obtained by solv- 
ing the Teukolsky radial equation and substituting 
a = 0, and they are clearly different from those in 
the Regge- Wheeler case. Hence, it is important to 
check whether both methods give the same results 
for QNM and with what precision. 

THE EPSILON METHOD 

For values of the parameters a, (3, 7, 8, r\ of general 
type, the confluent Heun function HcunC(a, f3, 7, 5, rj, z) 
has branching points in the complex z-plane at the sin- 
gular points z = 1 and z = 00. In the MAPLE package, as 
a branch cut is chosen the semi-infinite interval (1, 00) on 
the real axis. The presence of the branch cut may lead 
to the disappearance of some modes or their translation, 
since by changing the phase of the complex variable r, we 
may make a transition to another sheet of the multival- 
ued function. To avoid this, we use the epsilon method 
with which one can find the correct sheet and remain on 
it. This is done by introducing a small variation (| e |< 1) 



in the phase condition arg(r)+arg(aj) = — tt/2 (defined by 
the direction of steepest descent, see [IT]), with which 
one can move the branch cuts farther from the roots and 
thus avoid the jump discontinuity in the function. For 
more information on the epsilon method and the numer- 
ical procedures, see [20]. 

NUMERICAL RESULTS 

From the angular equation (1), it is clear that it can 
be solved explicitly without solving the system (1) and 
(2) and the values of I are known: i = 2,3,,.., In this 
paper, only the first value, I = 2, is used to find the 
QNMs with both radial equations. One can then either 
solve only the radial equations or solve the systems (1) 
and (2) with the appropriate values of the parameters. 
If one solves the problem as a two-dimensional system, 
making calculations with 15 digits of precision (and 32 
software floating-point digits), one obtains as expected, 
I = 1.99(9) + 1 x 10~ 17 i with the first digit different from 
digit 9 being the 17th. 

The numerical results for the frequencies are summed 
in Table I. 



n 


id from the Regge- Wheeler Eq. 


td from the Teukolsky Eq. 





0.7473433688+0. 1779246316i 


0.7473433676+0.1779246260i 


1 


0.6934219937+0. 5478297504i 


0.6934219698+0.5478298839i 


2 


0.6021069092+0.9565539668i 


0.6021069568+0.9565538786i 


3 


0.5030099245+1. 4102964056i 


0.5030097036+1.4102966442i 


4 


0.4150291600+1. 8936897821i 


0.4150291670+1.8936897747i 


5 


0.3385988052+2.3912161094i 


0.3385987682+2.3912160831i 


6 


0.2665046794+2.8958212549i 


0.2665047149+2.8958212406i 


7 


0. 1856446653+3.4076823515i 


0. 1856446394+3.4076823843i 


8 


-0.0306490371+3.9968237195i 


-0.0306490242+3. 9968236554i 


9 


0. 1265269702+4.6052896060i 


0.1265270059+4.6052895329i 


10 


0.15310679658+5. 1216534769i 


0.1531069231+5. 121653227H 



Table I. A list of the frequencies obtained for the QNMs of 
Schwarzschild black hole using the Regge- Wheeler equation 
and the Teukolsky equation. The modes with n < 5 are found 
for e = 0, modes from n > 5 - with e = —0.3. The first 5 
frequencies (n = — 4) were obtained also by Fiziev in [17] 
using exact solutions of RWE in the Heun functions 

From the table, one can see that the frequencies from 
the two types of equations coincide with at least 6 digits. 
A comparison between the RWE frequencies and the ones 
published by Andersson [21], published in [20] shows that 
the difference between the two results is smaller than 
5x 10~ 8 in most cases and is due to the numerical reasons. 

There are two important results from this study. First, 
as seen from Table I for both the RWE and the TRE, 
the mode number 8 has a small but nonzero real part. 
According to Leaver's evaluations this mode should be 
equal to + 3.998000i [15], with an exactly zero real 
part, if it is to correspond to the so-called algebraically 
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special mode. 

Algebraically special (AS) modes have a special place 
in the QNM studies [1]. The Andersson method is not 
applicable for them and these are excluded from his con- 
sideration. Berti, Cardoso and Starincts ([4, 5]) make 
a review on the results so far concerning these modes. 
Theoretically the 9th mode (n = 8) should be purely 
imaginary with value Ai, if it indeed corresponds to the 
AS case. In our results, even though purely imaginary 
modes do not pose a problem for the method, the real 
part of the 9th mode is distinctly not zero, and it has 
at least 7 stable digits when changing e in the interval 
discussed below for both RWE and TRE. This clearly 
shows that this mode does not agree with the hypothesis 
for the AS mode, which is to be expected since the AS 
mode should correspond to different boundary conditions 
- those of the so-called totally-transmission modes ([22]. 

The second important result is the dependence of the 
frequencies uj n on e. The direction of steepest descent is 
supposed to be the optimal direction in which the solu- 
tions satisfy the black hole boundary conditions on infin- 
ity in the first term approximation for asymptotic series 
for the Heun functions [17]. The validity of steepest de- 
scent method in its simplest form for the radial equations 
(2) in both cases under variations in this condition, how- 
ever, is still an open problem studied here for the first 
time. 

Using the e method, one can explore the intervals for 
e in which each mode can be found. The results for both 
RWE and TRE, as expected, coincide. Generally, the in- 
tervals into which each mode can be found narrow down 
when increasing n. While for the first 5 modes it is pos- 
sible to find u>„ = ±|3?(cj„)| + Q(uj n )i for positive and 
negative values of e in a certain interval , for n > 4, 
(but n / 8) the modes with a positive real part can be 
found only for negative values of e, and the dependency 
becomes w„(e) = — sgn(e)|5R(a; n )| + ^s(u! n )i. 

For n = 8, the mode has different behavior with 
respect to e - for e G [—0.75,-0.1], one finds a 
mode with negative real part and vice versa: (uj n= s = 
sgn(e) 0.030649006 + 3.996823690i). 

The so- found relation w„(e) needs to be examined fur- 
ther. For the case n = 8, similar (to some extent) be- 
havior was mentioned also in [22, 23] (and discussed in 
[5]). It was suggested that there are two AS modes which 
are symmetrical to the imaginary axis and perhaps may 
be related with the branch cut in the asymptotic of the 
RWE potential when u> is purely imaginary. Using the 



1 The ranges where each mode is found depend on e as follows: for 
n = 0, e e [=F0.8,±0.75], for n = 1, e 6 [t-0.8, ±0.45], for n = 
2, e G [=F0.8,±0.25], for n = 3, e e [^0.8, ±0.1], for n = 4, e € 
[=F0.8,0], where the first sign corresponds to frequencies with a 
positive real part and the second sign to those with negative real 
parts. The imaginary parts for each mode n coincide. 



e method applied on the asymptotics of the confluent 
Heun functions, one can directly obtain the place of the 
branch cut on the real axis as a function of e and they 
can be easily visualized plotting the solution . There- 
fore, the use of the confluent Heun functions and the e 
method offers a direct way to examine the solutions and 
their properties in relation to the branch cut in the com- 
plex r-plane, something that cannot be readily done in 
the continued fractions method generally used to obtain 
the QNMs. 

Further exploration of the dependence 5R(w n )(e) (or 
5(w n )(e)) in the intervals mentioned above shows that, 
for both the RWE and the TRE, it is approximately a 
periodic function with amplitude A and period L which 
change with n in a nontrivial way (Fig.l and Fig. 2). 
For n < 4, from the RWE and the TRE one obtains 
A T re ~ 10~ 6 « 10 3 x Ar W e, L rwe « L T re ~ 0.4 
and those values remain approximately constant with re- 
spect to n (n < 4). For n > 4, the dependence of A and 
L on n becomes more pronounced: the amplitudes and 
the periods of the RWE increase with n until they reach 
Arwe ~ 10" 6 , L WRE « 0.6 for n = 10. For the TRE, 
the amplitude and the period decrease to Atre ~ 10~ 8 , 
Ltre ~ 0.05. For n = 8, the two periodic behaviors have 
approximately equal amplitudes w 10 -7 . Those results 
hint that, although the so-obtained frequencies are stable 
with at least 6 digits with respect to e, there is also some 
finer dependence, the origin of which should be carefully 
investigated. 



CONCLUSION 

In this paper, were presented the QNMs for a 
Schwarzschild BH obtained from the RWE and the TRE, 
by solving the differential equations analytically in terms 
of confluent Heun functions. The QNMs from the TRE 
for the case s = — 2 were calculated for the first time and 
were found to coincide with the well-known QNMs from 
the RWE with precision of 6 digits. 

We demonstrated a new method for studying the sta- 
bility of the QNM calculations. The results show nontriv- 
ial dependence on small variation in the phase condition 
(the e method) which requires additional investigation. 

For the first time, the mode n = 8 was obtained di- 
rectly from the spectral condition on the exact analytical 
solutions of RWE and TRE and was found to have a 
nonzero real part, which proves that this mode is not the 
algebraically special mode. The mode in question is sta- 
ble with 6 digits of significance with respect to changes 
in e, which proves that its real part is indeed not zero. 

Those results presented here show the strength of us- 
ing confluent Heun functions to find QNMs of nonrotat- 
ing BHs and are encouraging in continuing this work in 
finding QNMs of rotating BHs. 
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Figure 1. Complex plots of the scaled QNMs from the two 
equations in the appropriate intervals for e: a)RWE n = 0..3 
b)TRE n = 0..3 c) RWE n = 4.. 10 d)TRE n = 4.. 10. Clearly, 
while for n — 0..3 the QNMs from the two equations give sim- 
ilar results, for n > 4, the variations in the frequencies from 
TRE happen on a much smaller scale and appear chaotic. 
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Figure 2. Plots of abs(uj)(e), where the solid line corresponds 
to the TRE modes and the dots - to the RWE modes for 
a) n=3, b) n=8, c) n=10. One can see the evolution of the 
amplitudes and the periods in both cases, when n increases 
from 3 to 10. 
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